function [error,alpha,beta,constant]=fn_modelfit_5d_PlotPoints(patient, slice,j,k,fig_num);
    
    %------------------------------------------------------------------------
    %   This file is part of the
    %   5D-Novel4DCT Toolbox  ("Novel4DCT-Toolbox")
    %   DH Thomas, Ph.D
    %   University of California, Los Angeles
    %   Contact: mailto:dhthomas@mednet.ucla.edu
    %------------------------------------------------------------------------
    % $Author: DHThomas $	$Date: 14-Apr-2014 15-23-24 $	$Revision: 0.1 $
    
    set(0,'defaultlinelinewidth',2)
    set(0,'defaultaxeslinewidth',2)
    
    if exist([patient.model_folder '/simple_coefs.mat'])>0;
        load([patient.model_folder '/simple_coefs.mat']);
        patient.volume_calibration  = simple_coefs;clear simple_coefs
        plot_calibrated = 1;
    else
        patient.volume_calibration  = [1,0]
        plot_calibrated = 0;
    end
    
    vox = sub2ind([patient.dim(2),patient.dim(3)],round(j),round(k));
    %   matlabpool 8
    %  h_temp=load('h_BW');
    %
    % h=h_temp.h_BW;clear h_temp;
    figure(fig_num);
    if length(findall(gcf,'type','axes'))>1;
        cla(subplot(2,2,2));
        cla(subplot(2,2,3));
        cla(subplot(2,2,4));
    end
    
    set(gcf,'units','normalized','position',patient.default_figure_position);
    h1 = subplot(2,2,1)
    imagesc(flipud(rot90(squeeze(patient.static(slice,:,:)))), [-950 350])
    set(gca,'YDir','normal')
    
    colormap gray
    axis image
    drawnow
    
    hold on;
    plot(j,k,'ro')
    
    
    cp=1;
    msx=15;
    mso=10
    
    %     indexs_vec=repmat(kk,size(jj,2),1);
    %     indexs_vec=indexs_vec(:)';
    
    
    
    [~,grid_z]=meshgrid(1:patient.dim(2),1:patient.dim(3));
    grid_z=permute(repmat(grid_z,[1,1,patient.scans]),[3,2,1]);
    
    
    [ dvf]=read_dvf_elastix_slice_toolbox(patient.dvf_folder,slice);
    
    voxel_w=double(bsxfun(@plus,grid_z,dvf.w(:,1:patient.dim(2),1:patient.dim(3))));
    voxel_w=voxel_w(:,:);
    
    for scan=patient.run_scans;
        %      hu_volt_tmp_s= ScaleTime(bellows_volt(kk,scan)',voxel_w(scan,:)');
        X1(scan,:) = nakeinterp1([1:patient.dim(3)]',patient.bellows_volt_drifted(:,scan),voxel_w(scan,:)'); %%% DOESN'T WORK - FIX THIS;
        
        X2(scan,:) = nakeinterp1([1:patient.dim(3)]',patient.flow_drifted(:,scan),voxel_w(scan,:)');
        %         X3(scan,:) = nakeinterp1(kk',ecg_index(kk,scan),voxel_w(scan,:)');
    end
    
%     X2(16,vox) = -X2(16,vox) ;

    static_mask_slice=squeeze(patient.static_mask(slice,:,:));
    static_mask_slice=static_mask_slice(:);
    
    dvf_u_vec=single(dvf.u(:,:));
    % dvf_u_vec=dvf_u_vec(:,:);%dvf_u_vec(leave_out,:)=[];
    dvf_v_vec=single(dvf.v(:,:));
    % dvf_v_vec=dvf_v_vec(:,:);%dvf_v_vec(leave_out,:)=[];
    dvf_w_vec=single(dvf.w(:,:));
    % dvf_w_vec=dvf_w_vec(:,:);%dvf_w_vec(leave_out,:)=[];
    
%     alpha  = load([patient.model_params_folder sprintf('/alpha_%d',slice)]);alpha = permute(alpha.data,[3,1,2]);alpha = alpha(:,:);
%     beta  = load([patient.model_params_folder sprintf('/beta_%d',slice)]);beta = permute(beta.data,[3,1,2]);beta = beta(:,:);
%     constant  = load([patient.model_params_folder sprintf('/constant_%d',slice)]);constant = permute(constant.data,[3,1,2]);constant = constant(:,:);
%     
%     tempfit_variables_x(:,:)=[constant(1,:); alpha(1,:) ;beta(1,:)];%E_vox\dvf_u_vec(patient.run_scans,vox);
%     tempfit_variables_y(:,:)=[constant(2,:); alpha(2,:); beta(2,:)];%E_vox\dvf_v_vec(patient.run_scans,vox);
%     tempfit_variables_z(:,:)=[constant(3,:) ;alpha(3,:); beta(3,:)];%E_vox\dvf_w_vec(patient.run_scans,vox);

    
    constant=zeros(patient.dim(2),patient.dim(3),3);
    alpha=zeros(patient.dim(2),patient.dim(3),3);
    beta=zeros(patient.dim(2),patient.dim(3),3);
    % gamma=zeros(length(jj),length(kk),3);
    
    cc=[1,0,0;0,1,0;0,0,1];
    cp=1;
    % for vox=1:size(ii,2)*size(jj,2)*size(kk,2);%size(X1,2);
    if static_mask_slice(vox)==0;
        display('WARNING: Voxel is outside of lung')
    end
    
    E_vox=[ones(size(X1(patient.run_scans,:),1),1) X1(patient.run_scans,vox) X2(patient.run_scans,vox)];
    
    tempfit_variables_x(:,vox)=E_vox\dvf_u_vec(patient.run_scans,vox);
    tempfit_variables_y(:,vox)=E_vox\dvf_v_vec(patient.run_scans,vox);
    tempfit_variables_z(:,vox)=E_vox\dvf_w_vec(patient.run_scans,vox);
    

    yfit_x(:,vox)=tempfit_variables_x(:,vox)'*E_vox';
    yfit_y(:,vox)=tempfit_variables_y(:,vox)'*E_vox';
    yfit_z(:,vox)=tempfit_variables_z(:,vox)'*E_vox';
    
%     yfit_x2(:,vox)=tempfit_variables_x2(:,vox)'*E_vox';
%     yfit_y2(:,vox)=tempfit_variables_y2(:,vox)'*E_vox';
%     yfit_z2(:,vox)=tempfit_variables_z2(:,vox)'*E_vox';

    error_u_vox=yfit_x(:,vox)-dvf_u_vec(patient.run_scans,vox);
    error_v_vox=yfit_y(:,vox)-dvf_v_vec(patient.run_scans,vox);
    error_w_vox=yfit_z(:,vox)-dvf_w_vec(patient.run_scans,vox);
    error_vec_vox_scans=bsxfun(@hypot,error_u_vox,bsxfun(@hypot,error_v_vox,error_w_vox));
    
    error_vec_vox(cp)=mean(bsxfun(@hypot,error_u_vox,bsxfun(@hypot,error_v_vox,error_w_vox)));
    

    
    figure(fig_num);
    h2 = subplot(2,2,2)
    hold on
    p1=plot(dvf_v_vec(patient.run_scans,vox),dvf_w_vec(patient.run_scans,vox),'x','MarkerSize',msx);
    p12=plot(dvf_v_vec(patient.ref,vox),dvf_w_vec(patient.ref,vox),'kx','MarkerSize',msx);
    
    p2=plot(yfit_y(:,vox),yfit_z(:,vox),'ro','Color',cc(cp,:),'MarkerSize',mso);
    p22=plot(yfit_y(patient.ref,vox),yfit_z(patient.ref,vox),'ko','MarkerSize',mso);
    
    plot([dvf_v_vec(patient.run_scans,vox),yfit_y(:,vox)]',[dvf_w_vec(patient.run_scans,vox),yfit_z(:,vox)]','k--','Color',cc(cp,:));
    hold off
    ylabel('Craniocaudal [mm]')
    xlabel('Anteposterior [mm]')
    %         axis([-1.5 4 -4 3])
    axis image
    title(sprintf('Mean Error = %.2f, Max Error = %.2f [mm]', error_vec_vox, max(error_vec_vox_scans)))
    %     legend([p1(1,1);p2(1,1)],'Measured','Fit')
    N=30
    set(gca,'FontSize',15)
    h_xlabel = get(gca,'XLabel');
    set(h_xlabel,'FontSize',N);
    h_ylabel = get(gca,'YLabel');
    set(h_ylabel,'FontSize',N);
    
    h3 = subplot(2,2,3)
    hold on;grid on
    p1=plot3(polyval(patient.volume_calibration, X1(patient.run_scans,vox)),  X2(patient.run_scans,vox),dvf_w_vec(patient.run_scans,vox),'x','MarkerSize',msx);
    %         p12=plot(polyval(patient.volume_calibration, X1(patient.patient.ref,vox)),dvf_w_vec(patient.patient.ref,vox),'kx','MarkerSize',msx);
    p2=plot3(polyval(patient.volume_calibration, X1(patient.run_scans,vox)), X2(patient.run_scans,vox),yfit_z(:,vox),'ro','Color',cc(cp,:),'MarkerSize',mso);
    plot3([polyval(patient.volume_calibration, X1(patient.run_scans,vox)),polyval(patient.volume_calibration, X1(patient.run_scans,vox))]',[X2(patient.run_scans,vox),X2(patient.run_scans,vox)]',[dvf_w_vec(patient.run_scans,vox),yfit_z(:,vox)]','k--','Color',cc(cp,:));
    hold off
    if plot_calibrated>0
        xlabel('Volume [L]')
        ylabel('Flow [L/s]');
    else
        xlabel('Bellows Voltage [V]')
        ylabel('Flow [V/s]');
    end
    
    zlabel('Craniocaudal [mm]')
    %          legend([p1(1,1);p2(1,1)],'Registration','Model')
    
    
    set(gca,'FontSize',15)
    h_xlabel = get(gca,'XLabel');
    set(h_xlabel,'FontSize',N);
    h_ylabel = get(gca,'YLabel');
    set(h_ylabel,'FontSize',N);
    h_zlabel = get(gca,'ZLabel');
    set(h_zlabel,'FontSize',N);
    view(3)
    
    
    h4 = subplot(2,2,4)
    
    hold on
    p1=plot(polyval(patient.volume_calibration, X1(patient.run_scans,vox)),dvf_w_vec(patient.run_scans,vox),'x','MarkerSize',msx);
    p12=plot(polyval(patient.volume_calibration, X1(patient.ref,vox)),dvf_w_vec(patient.ref,vox),'kx','MarkerSize',msx);
    
    p2=plot(polyval(patient.volume_calibration, X1(patient.run_scans,vox)),yfit_z(:,vox),'ro','Color',cc(cp,:),'MarkerSize',mso);
    plot([polyval(patient.volume_calibration, X1(patient.run_scans,vox)),polyval(patient.volume_calibration, X1(patient.run_scans,vox))]',[dvf_w_vec(patient.run_scans,vox),yfit_z(:,vox)]','k--','Color',cc(cp,:));
    hold off
    if plot_calibrated>0
        xlabel('Volume [L]')
    else
        xlabel('Bellows Voltage [V]')
    end
    
    
    ylabel('Craniocaudal [mm]')
    legend([p1(1,1);p2(1,1);p12(1,1)],'Registration','Model','Reference')
    
    set(gca,'FontSize',15)
    h_xlabel = get(gca,'XLabel');
    set(h_xlabel,'FontSize',N);
    h_ylabel = get(gca,'YLabel');
    set(h_ylabel,'FontSize',N);
    

    
    error = nanmean(error_vec_vox_scans);
    
    
    